% m file generated from linear noise model SIR
function [T,Y] = prova(tf,points)
if(nargin==0)
	tf=100; points=1000;
elseif (nargin == 1)
	points = 1000;
end;

%model parameters
par(1) = 1000.0;  % N
par(2) = 0.99;  % s0
par(3) = 1.0;  % normalise
par(4) = 10.0;  % ki
par(5) = ( ( par(3) == 1 ) * par(4) + (1 - ( par(3) == 1 ) ) * ( par(4) / par(1) ) );  % kiN
par(6) = 0.1;  % kr
par(7) = 0.25;  % ks


%model variables
var(1) = ( ( par(3) == 1 ) * par(2) + (1 - ( par(3) == 1 ) ) * ( par(2) * par(1) ) );  % S
var(2) = ( ( par(3) == 1 ) * ( 1 - var(1) ) + (1 - ( par(3) == 1 ) ) * ( par(1) - var(1) ) );  % I
var(3) = 0.0;  % R
var(4) = 0.0;  % var.S
var(5) = 0.0;  % cov.S.I
var(6) = 0.0;  % cov.S.R
var(7) = 0.0;  % var.I
var(8) = 0.0;  % cov.I.R
var(9) = 0.0;  % var.R


%defined function


%ode function
function dx = odefun(t,var)

n = 3;
m = 3;
code = @(i,j)(n*i - (i-1)*(i-2)/2 + j - i + 1);

f = [ 	( ( par(5) * var(1) ) * var(2) );
	( par(6) * var(2) );
	( par(7) * var(3) )];

S = [ 	-1.0, 0, 1.0;
	1.0, -1.0, 0;
	0, 1.0, -1.0 ];

Jf = zeros(m,n);
Jf(1,1) = ( par(5) * var(2) );
Jf(1,2) = ( par(5) * var(1) );
Jf(1,3) = 0.0;
Jf(2,1) = 0.0;
Jf(2,2) = par(6);
Jf(2,3) = 0.0;
Jf(3,1) = 0.0;
Jf(3,2) = 0.0;
Jf(3,3) = par(7);

J = S*Jf;

C = zeros(3);
for i=1:n
	for j=i:n
		C(i,j) = var(code(i,j));
		if (i ~= j)
		C(j,i) = var(code(i,j));
		end;
	end;
end;
dx = S*f;
dC = J*C + C*J' + S*diag(f)*S';
for i=1:n
	for j=i:n
	dx(code(i,j)) = dC(i,j);
	end;
end;
end


%solving ode
delta = tf/points;
tspan = 0:delta:tf;
[T,Y] = ode45(@odefun,tspan,var);

end

